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Analysis  of  Smoke  Trail  Photographs  to  Determine 

Stratospheric  Winds  and  Shears 

1,  INTRODUCTION 

Detailed  measurements  of  the  vertical  wind  profile  at  stratospheric  altitudes 
are  required  to  evaluate  atmospheric  transport  processes  of  importance  in  estab- 
lishing dispersal  rates  of  SST  emissions.  These  rates  are  controlled  by  the  inten- 
sity of  turbulence  and  the  local  shear  fields  at  spatial  scales  of  the  order  of  a few 
meters.  Conventional  balloon  observations  and  meteorological  rockets  are  not 
directed  towards  the  definition  of  flow  fields  in  the  atmosphere  at  such  small  scales. 
To  obtain  the  necessary  high  resolution,  we  have  relied  on  observations  and  analysis 
of  visible  tracer  trails  deposited  from  rockets  programmed  to  release  a dense 
TiCl^  smoke  from  the  time  the  rocket  reaches  an  altitude  of  some  15  km  to  apogee. 
The  trails  are  photographed  at  regular  intervals  from  at  least  two  different  locations 
to  provide  positional  and  timing  information  for  triangulation,  after  which  the  wind 
field  parameters  and  the  shears  are  computed. 

2.  OPERATIONAL  CONCEPT 

To  obtain  the  wind  field  parameters  photogrammetric  measurements  are  needed 
on  photographs  of  successive  positions  of  the  trails  laid  by  the  rocket  vehicle.  For 
each  field  program,  photographs  were  taken  from  three  different  sites  at 

(Received  for  publication  8 October  1976) 
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White  Sands  Missile  Range.  During  February  1973,  for  example,  measurements 
were  made  by  triangulation  cameras  located  at  Dona  Ana,  "A"  Mountain  and 
Two  Buttes.  For  the  June  1973  and  June  1975  programs,  the  sites  selected  were 
Two  Buttes,  Hotel  and  Seehorn.  Figure  1 shows  their  locations  as  well  at  the 
location  of  the  launching  complexes.  Coordinates  of  these  stations  are  listed  in 
Appendix  A.  At  two  of  these  sites  three  cameras  were  used.  The  third  site  had 
two  cameras.  In  all  cases  a camera  equipped  w>th  either  a 20-in.  or  7-in.  focal 
length  lens  was  co-mounted  with  a similar  camera  having  a 7-'n.  lens.  The  field 
of  view  of  the  20-in.  lens  was  12  X 12  deg  and  was  adjusted  to  lie  in  the  upper  third 
of  the  field  of  the  7-in.  camera  which  viewed  a 36  X 36  deg  field.  At  the  sites  where 
a third  camera  was  used  it  carried  a 7-in.  lens  and  it  was  mounted  on  a separate 
rig  oriented  to  photograph  essentially  the  same  field  seen  by  the  7 -in.  co-mounted 
camera.  The  operation  of  the  camera  and  film  processing  details  are  given  in 
Appendix  B.  Figure  2 shows  the  evolution  of  a trail.  The  gaps  are  created  inten- 
tionally to  provide  easily  identifiable  end  points  on  views  of  the  trail  photographed 
from  different  sites.  This  procedure,  introduced  in  1975,  facilitates  the  triangula- 
tion particularly  after  fine  structure  appears  and  the  trail  image  has  developed  self 
intersections.  As  can  be  seen  from  Figure  2 the  trail  image  becomes  rather  tangled 
in  a surprisingly  short  time.  Each  of  the  triangulation  photographs  measures  12.  5 
X 12.5  cm  and  contains  the  trail  image  plus  the  images  of  fiducial  marks  that, 
together  with  the  normal  to  the  film  plane,  define  a local  Cartesian  system  of  coor- 
dinates. The  orientation  of  this  local  system  with  respect  to  some  reference  sys- 
tem (for  example,  horizon,  geodetic,  etc.)  is  determined  by  means  of  auxiliary 
photographs  of  the  sky  taken  at  known  times  during  the  night.  Measurements  of 
star  image  coordinates,  referred  to  the  film  plane  axes,  are  used  for  tnis  purpose. 

If  the  cameras  are  not  moved  from  the  time  the  star  background  is  photographed, 

a single  transformation  can  be  found  to  change  the  film  plane  coordinates  into  angu- 

1 2 

lar  coordinates  to  identify  the  orientation  of  lines  of  sight  to  points  on  the  trails.  ’ 

If  the  cameras  are  moved,  the  coordinates  of  fixed  objects  that  appear  on  photo- 
graphs taken  before  and  after  the  camera  realighment  (such  as  mountain  tops  or 
specially  constructed  targets)  are  needed  to  define  additional  transformation  of 
coordinates.  Computations  have  shown  that  in  order  to  obtain  a spatial  reso- 
lution of  the  order  of  10  m at  a distance  of  some  30  km  it  is  necessary  to  sub- 
ject the  trail  photographs  to  a scanning  operation  on  densitometric  equipment 
capable  of  generating  a raster  with  consecutive  line  separation  of  the  order  of 
50  n used  in  combination  with  scanning  apertures  approximately  50  square. 


1.  Quesada,  A.  F.  ( 1 97  1)  AFCRL-7  1 -0233,  ERP  No.  351. 

2.  Quesada,  A.  F.  (1975)  AFCRL-TR-75-045 1 , ERP  No.  527. 
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Figure  1.  Location  of  Triangulation  Sites 


Experience  has  shown  that  many  trail  sections  can  be  scanned  at  half  this  resolu- 
tion. Occasionally,  there  are  film  areas  where  the  appearance  of  very  fine  detail 
requires  a scan  raster  and  aperture  of  25  microns.  In  any  event,  the  positional 
accuracy  of  the  system  must  be  such  that  x and  coordinates  of  points  along  the  trail 
are  read  with  errors  not  exceeding  one  to  two  fx /centimeter.  Appendix  C describes 
a pictorial  analysis  facility  (PAF)  that  meets  these  requirements.  In  order  to 
preserve  the  high  accuracy  derived  from  scanning  the  trail  photographs  with  the 
PAF  equipment,  the  positional  information  is  digitized  and  recorded  on  magnetic 
tape.  The  data  on  the  tape,  consisting  of  star  coordinates,  fiducial  mark  coor- 
dinates, trail  coordinates  and  (possibly)  auxiliary  target  coordinates,  are  processed 
on  the  CDC  6600  Computer  at  AFGL  to  determine  the  successive  positions  of  points 
on  the  trails  at  altitude  increments  of  10  m.  From  three  site  positional  data,  wind 
velocities  and  error  estimates  are  computed  as  functions  of  altitude  and  time. 


3.  TRI  ANGULATION  PROCEDURE 

The  assignment  of  space  coordinator  to  each  point  on  the  trail  is  a function  of 
the  orientation  of  the  cameras  and  of  the  film-plane  coordinates  of  each  point.  When 
the  star  calibration  data  are  reduced,  it  is  found  that  camera  orientation  param- 
eters and  camera  focal  length  can  be  determined  with  errors  measured  by  standard 
deviations  of  0.  01  to  0.  03  deg  and  0.  01  cm,  respectively. 

The  triangulation  itself  uses  as  a criterion  for  matching  points  on  photographs 
of  the  trail  corresponding  to  different  sites  the  equality  of  dihedral  angles  formed 
by  a line  connecting  a pair  of  observation  sites  and  the  line  of  sight  vectors  to  a point 
on  the  trail.  The  triangulation  package  consists  of  a set  of  programs  to  implement 
the  following  procedures: 

(1)  Computation  of  camera  focal  length  using  star  positions  on  the  film  plane, 

(2)  Determination  of  camera  orientation  from  star  positions, 

(3)  Determination  of  camera  orientation  from  coordinates  of  auxiliary  targets, 

(4)  Correction  for  atmospheric  refraction  when  elevation  angles  are  below 
10  deg, 

(5)  Determination  of  coordinates  of  points  along  the  trail  from  two-site 
positional  data, 

(6)  Determination  of  statistics  of  point  coordinates  along  the  trail  when  three- 
site  data  are  available, 

(7)  Conversion  to  coordinate  systems  more  suitable  to  specify  location  of  trails 
in  space, 

(8)  Computation,  tabulation  and  graphical  displays  of  winds,  wind  shears  and 
their  errors. 

Some  of  these  programs  are  listed  in  Appendix  D. 
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Figures  3 and  4 respectively,  show  the  variation  in  the  E-W  and  N-S  compo- 
nents of  the  wind  for  one  of  the  trails  released  during  the  June  1973  experiments, 
computed  from  two-site  triangulation  data. 


-5.0  0.0  5.0  10.0  15-0  ;0.0  JS.O  3C-C  35.0 
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Figure  3.  East-West  Component  of  Wind  Velocity 
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Figure  4.  North-South  Component 
of  Wind  Velocity 
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The  early  triangulation  attempts  and  the  subsequent  calculation  of  winds  and 
shears  did  not  lead  to  satisfactory  results.  Lack  of  consistency  was  found  in  three 
areas: 

(1)  The  position  in  space  of  a given  point  along  the  trail,  obtained  from  data 
corresponding  to  a pair  of  sites,  did  not  match  the  position  of  the  same  point  when 
it  was  calculated  from  the  data  corresponding  to  a different  pairing  of  sites.  The 
ground  plane  discrepancies  for  some  points  were  of  the  order  of  50  to  200  meters. 

(2)  The  wind  profiles  corresponding  to  different  pairings  did  not  agree,  the 
discrepancies  in  this  case  being  small  shifts  in  the  altitudes  assigned  to  prominent 
features  of  the  profile,  and  marked  differences  in  the  fine  structure  of  the  profiles. 

(3)  The  shears  showed  abnormally  high  values  at  too  many  points  in  the  altitude 
range  that  was  covered. 

An  intensive  search  for  possible  causes  for  these  discrepancies  led  to  the  con- 
clusion that  small  changes  in  the  orientation  parameters  of  the  camera  and  in  the 
camera  focal  length  were  at  fault.  These  changes  are  for  the  most  part  produced 
by  agents  not  under  the  control  of  the  operator.  It  was  found,  for  example,  that 
camera  mounts  can  move  slightly  when  the  film  magazines  are  taken  on  and  off  to 
protect  the  film  against  excessive  heating  during  the  day.  They  can  also  move  when 
the  canvas  covers  are  put  on  or  taken  off.  The  focal  length  changes  during  the  ex- 
periment due  to  heating  of  the  cameras  as  the  sun  rises,  etc.  All  of  these  factors 
can  produce  changes  whose  magnitude  is  a significant  fraction  of  a standard  devia- 
tion in  the  parameter  values,  that  is,  some  tens  of  p for  the  focal  length  and  a few 
tenths  of  a mrad  in  the  orientation  parameters. 

To  apply  corrections  to  the  computed  parameter  values  it  is  necessary  to  intro- 
duce an  optimization  procedure  which,  starting  from  the  triangulated  values  of 
positions  along  the  trail  corresponding  to  a given  paring  of  sites,  computes  the  film 
plane  coordinates  of  the  trail  for  the  site  that  is  to  be  corrected  and  compares  it 
with  the  original  film  plane  coordinates  that  served  as  input  data  for  the  triangula- 
tion. 

The  program  stop  when  the  sum  of  the  absolute  values  of  the  differences 
reaches  a minimum.  It  ordinarily  takes  two  or  three  iterations  to  zero  in  on  a set 
of  parameter  values  that  bring  positions  computed  from  different  pairings  to  better 
than  10  m (1  part  in  5000)  and  makes  the  wind  profiles  coincide  in  all  of  the  gross 
features  and  most  of  the  fine  structure.  The  shears  are  less  pronounced  and  large 
values  occur  at  very  few  altitudes. 

Figure  5 shows  the  magnitude  of  the  wind  vs  altitude  computed  from  different 
site  pairings  for  the  June  1975  trail  of  Figures  3 and  4.  Figure  6 displays  the 
corresponding  shears. 
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Figure  5.  Magnitude  of  Wind  Velocity  Computed  from  Two  Different  Site  Pairs 


SMEAR  f SEC«R-2  I 

Figure  6.  Profile  of  Wind  Shear  Squared 

4.  CONCLUSIONS 

It  has  been  established  that  stratospheric  winds  with  a 10  m altitude  resolution 
can  be  determined  from  time-lapse  photographs  of  smoke  trails.  The  orientation 
of  the  cameras  is  critical  and  must  be  measured  with  a precision  that  requires 
optimizing  the  focal  length  and  the  orientation  parameters  of  the  camera.  W hen  this 


is  done,  the  triangulation  procedure  applied  to  3 or  more  sites  leads  to  results  that 
are  consistent  in  assigning  the  space  coordinates  of  points  on  the  trail  as  well  as 
to  agreement  of  wind  profiles  in  both,  significant  features  and  fine  structure.  Com- 
plementary optimization  and  triangulation  are,  therefore,  a technique  of  data  reduc- 
tion which  for  the  first  time  has  made  possible  the  computation  of  winds  and 

3 

shears  with  a resolution  three  times  higher  than  has  hitherto  been  reported. 


3.  Miller,  R.  W. 


Henry,  R. M. 


and  Rowe,  IVI.G.  (1965)  NASA  TN  D-2937. 


Appendix  A 


Coordinates  of  Triangulation  Sites 


Station 

Latitude 

Dona  Ana 

32° 

08  ft 

58. 2340 

"A"  Mountain 

32 

17 

17.  2409 

Horn 

32 

44 

24.  1298 

Hotel 

32 

43 

47. 2904 

Two  Buttes 

32 

42 

12.  8100 

Seehorn 

32 

45 

20.  747 

T-5 

32 

25 

19. 0256 

VVSMR 

(Bldg.  21620) 

32 

21 

39. 6309 

Holloman  AFB 
(Bid.  1251) 

32 

55 

27. 0990 

Geodetic  Coordinates 
WSTM  Coordinates 


Longitude 

Altitude 

106° 

30  ft 

08.  2121  in. 

4.  078.  44  ft 

106 

41 

45. 9135 

4,  754. 0 

106 

29 

29.  1790 

4, 425. 42 

106 

12 

50. 7995 

3,  992.  52 

106 

07 

35.  1100 

4,  555.  6 

106 

29 

12.  164 

4,  343.7 

106 

32 

59.  4567 

5,404. 8 

106 

24 

32.  4840 

3,  937. 6 

106 

04 

16. 6859 

4,  185.  2 
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Appendix  B 

Camera  Operation  and  Film  Development 


Cameras 

Modified  K-46  aerial  cameras  were  equipped  with  either  20-in.  or  7-in.  focal 
length  lenses.  Pairs  of  20-in.  and  7-in.  or  7-in.  and  7-in.  cameras  were  arranged 
on  a single  mount  for  each  site  and  rigidly  locked  to  move  as  one  unit  with  their 
optic  axes  offset  in  elevation  by  about  7 degrees.  At  two  of  the  sites  an  additional 
7 -in.  camera  was  mounted  separately.  Star  calibrations  were  made,  when  the 
weather  permitted,  the  night  before  the  start  of  the  experiments,  and  the  night  be- 
fore each  of  the  following  launch  sequences.  Exposure  times  of  15,  30,  60  and 
120  sec  were  used.  The  elevation  of  the  7-in.  cameras  was  chosen  to  allow  full 
coverage  of  the  trails  and,  in  addition,  show  ground  features  or  special  targets  near 
the  bottom  edge  of  the  film.  Daytime  photography  of  the  ground  features  was  re- 
quested following  the  star  calibrations.  It  is  thus  possible  to  correlate  ground 
features  and  star  calibrations  to  determine  the  camera  orientation  following  possible 
changes  in  the  camera  azimuths.  Azimuth  adjustment  during  the  experiments  may 
be  necessary  to  center  the  trail  images  within  the  field  of  view  of  the  cameras.  All 
cameras  were  provided  with  Polaroid  filters  mounted  on  top  of  Wratten  29  filters. 
The  Polaroids  were  rotated  shortly  before  launch  until  minimum  sky  brightness 
was  observed. 


For  the  July  1976  experiments,  a star  calibration  was  made  about  1 hr  before  the 
trail  was  released  in  order  to  minimize  camera  movements  due  to  film  magazine 
loading  and  unloading. 


Preceding  page  blank 
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Film  Development 

The  negatives  (Kodak  2402)  were  developed  to  have  the  following  character- 
istics: 

(1)  Gamma  in  the  range  1.2  to  1.5,  determined  from  measurements  of  step 
wedge  images  exposed  through  a Wratten  29  filter  and  a Polaroid  filter. 

(2)  Exposure  for  the  sky  background  leading  to  a density  on  the  negative 
approximately  0.  3D  above  fog  was  desired.  In  case  the  exposure  varied  from  the 
center  to  the  edge  of  the  field,  the  above  criterion  was  to  apply  to  the  background 
density  near  the  edges  of  the  film,  to  provide  reference  lines  for  coordinate  mea- 
surements in  case  the  fiducial  marks  were  lost. 

(3)  A developer  was  specified  (D-76)  that  " as  compatible  with  the  contrast 
requirements  stipulated  above  and  the  fine  grain  to  be  expected  from  a film  rated 
ASA  125. 
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Appendix  C 


Film  Measurements 


Film  measurements  are  carried  out  under  contract  by  Information  Design,  Inc. 
of  Bedford,  Massachusetts,  using  pictorial  analysis  equipment  consisting  of  a 
photo  scanner  interfaced  to  a general  purpose  digital  computer  wh:ch  also  controls 
interactive  display  and  recording  peripheral  units.  Characteristics  of  this  equip- 
ment are  listed  below. 

(1)  P-1000  computer  controlled  Photo  Scanner  with  capability  of  scanning  at 
any  of  the  following  X-Y  rasters: 

(a)  12.5  jU 

(b)  25  |i 

(c)  50  (Ll 

(d)  100  M 

(e)  200  4 

Any  of  the  following  apertures  can  be  used  in  combination  with  any  of  the  above 
rasters: 

(a)  12.  5 /Lt2 

(b)  25  fx2 

(c)  50  M2 

(d)  100  H2 

(e>  200  p2 

(f)  500  /x  by  25  /Lt  slit 

Operation  of  the  Photo  Scanner  is  completely  computer  controlled. 


X raster  size,  Y raster  size,  and  the  x and  y boundary  lines  of  the  image  area 
to  be  scanned  are  controlled  by  the  software.  Each  aperture  is  centered  at  the 
intersections  of  the  X-Y  lattice  raster. 

(2)  Alpha- 16  general  purpose  digital  computer  interfaced  to  the  P-1000. 

(3)  Proprietary  interactive  display  interfaced  to  the  computer. 

The  interactive  display  includes  a cathode  ray  tube,  keyboard  and  special 
digital  knob  controls.  Using  the  knob  controls,  variables  such  as  cursor  positions 
or  image  brightness  thresholds  can  be  input  to  the  computer  program  with  13-bit 
accuracy.  The  system  provides  the  means  for  rapid  and  accurate  manipulation  of 
image  data. 

(4)  Proprietary  system  control  console  that  allows  flexible  interactive  control 
and  editing  of  data. 

(5)  Magnetic  tape  800  b.  p.  i.  , 9 tracks,  program  for  conversion  to  7 track  at 
556  or  800  b.  p.  i. 

(6)  High  speed  paper  tape  reader  and  punch. 

(7)  General  purpose  computer  software  adapted  to  the  above  system  configura- 
tion. 

(8)  High  speed  line  printer. 

The  Pictorial  Analysis  Facility  (PAF)  has  been  programmed  to  detect  and 
digitize  the  smoke  trail  images  automatically.  However,  an  operator  can  provide 
guidance  via  the  display  in  order  to  initialize  an  image,  track  exceedingly  thin 
trails  and  track  self  intersections  of  the  trail  image.  Positional  x,y  accuracy  of 
the  densitometer  is  specified  by  the  manufacturer  as  1 H rms/cm. 

Each  frame  is  scanned  perpendicular  to  the  direction  of  smoke  trail  ascent. 
Thus,  scan  lines  generally  cross  the  smoke  trail  images  at  approximately  right 
angles.  The  array  of  optical  density  measurements  for  each  scan  line  is  analyzed 
automatically  by  computer  program  to  find  the  excursion  caused  by  the  smoke  trail 
image. 

The  operator  monitors  the  digitization  of  the  initial  frame  of  each  sequence. 
After  these  are  digitized,  the  computer  retains  an  abbreviated  model  of  the  n-1 
frame  smoke  trail  image.  This  model  consists  of  smoke  trail  image  x,  y-coordin- 
ates  spaced  at  1 mm  distances  along  the  vertical  dimension  of  the  trail  image. 

When  the  nth  frame  is  digitized,  frame-to-frame  correlation  is  used  to  direct 
the  scanner  to  the  appropriate  location  of  the  film.  The  search  for  the  trail  on  the 
nth  frame  is  centered  on  the  location  of  the  trail  on  the  n-1  frame. 

In  addition,  scan-line  to  line-scan  correlation  is  used  to  compare  x,y  coor- 
dinate locations  for  the  n-1  and  nth  scan  lines  in  any  frame.  Using  information 
both  from  the  previous  scan  line  on  the  same  frame  and  f)'om  the  same  scan  line 
on  the  previous  frame  provides  a good  set  of  expected  coordinates  for  the  current 
scan  line.  Initially  the  operator  may  define  a band  of  allowable  deviation  from  these 
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coordinates,  or  the  program  may  use  preset  values.  PAF  scans  each  smoke  trail 
automatically  until  an  out-of-tolerance  coordinate  pair  is  detected  and  then  sum- 
mons the  operator  to  interactively  guide  the  scanning  process.  If  necessary,  the 
operator  may  direct  the  PAF  to  repeat  the  last  few  scan  lines.  The  excursion  in 
optical  density  values  caused  by  the  smoke  trail  images  is  analyzed  by  computer 
to  determine  the  center  of  the  image  for  each  scan  line.  The  results  of  this  analy- 
sis are  smoke  trail  coordinates  ready  for  storage  on  magnetic  tape. 

Trail  image  thickness  are  calculated  each  time  an  x,  y coordinate  pair  is  de- 
tected. Trail  image  thickness  of  less  than  a preselected  threshold  causes  the  com- 
puter to  summon  the  operator  for  supervision  and  interactive  guidance. 

At  this  point,  the  computer  displays  an  image  of  the  area  surrounding  the 
anticipated  location  of  the  smoke  trail  image  or  an  area  designated  by  the  operator. 
The  operator  guides  the  digitization  by  viewing  a cross  section  of  the  optical  density 
profile  for  the  smoke  trail  image,  and  moving  a cursor  to  designate  the  position  of 
the  smoke  trail.  When  the  trail  widens  beyond  the  threshold,  the  PAF  notifies  the 
operator  that  it  has  returned  to  the  automatic  mode  of  operation. 

All  f'ducial  mark,  star,  auxiliary  target  and  smoke  trail  coordinates  as  re- 
quired are  stored  on  magnetic  tape  fully  compatible  with  the  CDC  6600  computer 
system  currently  in  use  at  AFGL. 


Appendix  D 

Computer  Programs 


Preceding  page  blank 
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1 PRrGRAM  TRIANG(INPUT  ,CUTFUT»TAPE5  = INPLIT  ,TAP£5= OUTPUT  , 

lTArE2r. TAPE  1= 150, TAPE£  = 15C»TAFcil  = l5J»TAPEi2  = l5J,TAPE3) 

r 

THIS  F-OGRAM  TRIANGULATES  A TRAIL  RELEASE,  ANO  CALCULATES  THE 
c r GEOCENTRIC  ANO  GEODETIC  COORDINATES. 

r TH:  LIMITS  ON  THIS  PROGRAM  ARE,  3 SITES,  AND  3031  POINTS  PER  ^ITi. 


REAL  TOENT 

DIMENSION  XXC 3000) ,VY( 3000) 

10  COMMON  C ( 3,30  00) , IOENT  (7,4)  ,N < 31  ,RHO(3 , 3) , 0 ( 3)  , 

1AV2C3) ,OUC3),FL<3>  ,AZIM(3>, AtL  (3) ,DELTAE<3) , 

2TITLS  (3,3»  ,XI(9)  , YIC,)  , ID(9)  , 

3TEMP1 (3) ,TEMP2(3> • TEMP3 131 ,CC (3) , AV (3> 

9, Gh(3 ) , ACEK (3  > , DEC (3) 

15  PI=3. 1L159265353979 

RA0=18C./PI 
J=1 

REWIND  20 
REWIND  1 
TT  REWIND  2 

REWIND  3 
READ(5,350>  NT  1 
350  FORMAT (131 
WRITE (E , 1 0 ) 

25  10  FO  °M  A T (1H1) 

DO  300  Nr=l,NTl 
RFNINO  11 
REWIND  12 
J=1 

30  Ec  IUMT  = 10  + J 

D 

C READ  THE  SITE  NAME,  AND  ITS  GEOOETIC  COORDINATES.  (CAT, LON) 

C 

READ<5,30>  (IOENT (I  , J),I  = 1,2>  , I DENT  (3,  J)  , IDENT  (•»  ,J)  , I DENT  (7  ,J) 

35  30  FORMAT(2A10,2F1J.9,F10 .5) 

IF  ( IOENT (1,J) .EQ. 1 1 H ) GO  TO  200 

IOENT (3, J»=I0ENT(3, JI/RAO 
ID  "NT  (-  , J ) = IOENT  ( 9 , J ) / RAD 
TEMPI (1) =IO ENT (3,  J)*RAO 
90  TEMPl(2)=IOENT (9,J)*RA9 

PRINT  50,  (IDENT(I,J> ,1  = 1,2) , TE MP1 ( 1 ) , TEMPI < 2 > , IDENT C 7, J) 

50  FORMAT (//  ,58X,2A1C,F10.L,*N.  LAT. *F10 . 4*W.  LONG.*/ 

1 56X,F7.2*METERS*) 

D ' ' 

95  C CONVERT  GEOOETIC  TO  GEOCENTRIC  COORDINATES, 

r; 

CALL  CONVRT (1 , IOEN  T ( 3 , J»  , 1DEN  T (9  , J > , IOENT (7 , J> /10  0 C . , IOENT (5 , J ) , 
1 IOENT (6 , J) , RHO ( 1 . J ) ) 

C 

50  c.  seao  The  universal  Time  and  date. 

r 

RE  AO (5,1 20)  JHR,JMIN,SEC,MO,IDAY,IYR 
120  FOPMAT  (2I3,F7.0,7X,I2,2I3) 

U T = ( JHP ♦ (JNIN*SEC/60.)/60.) /29. 

55  C 

C CALCULATE  SIOEPIAL  TIME. 


j 


IP| |g 


I 


I 


C 3 LL  JULOAY (JDAY, IYR,Mj,IOAY) 

F 0 TU=(JCAY*3. 5+UT-2 415020. >/36525. 

PIJ  = lt-  .♦(  38.+45.S36/6  0,  ) / 60.  ♦ (d64013  4.54  2*T'J+  0.  -529*TU*T'J)  / 360  . 

=?UU=AM"-0  <Ril  ,2-» . > 

GM4:(  13.  »UT*24.*RUU)  / 2-*,*2.  0*PI 
IP (GHA. G:.2.0*PI)  GHA=GHA-2.0*PI 

c r « 

r RrAO  TMI  camera  focal  length  in  cm. 

r 

REAO<5,20)  FLU) 

W3ITC(  ,110)  FL(J) 

T5  113  FORMAT  C r X , *FOCAL  length  =*FT  , 3*  CM*,/) 


P E AO  CJNGTANTS  WHICH  OEFINt  THE  CAMERA  OR IEN T A TIO N . < IN  OEGRuEG) 

r 

RE  AO ( r , 2 0 ) AZIMM,EL£V,CELCAE< J) 

75  PRINT  360 

56C  FORMAT (IX,*  AZ  EL  TILT*) 

PRINT  P0,AZIMM,ELEY,D£LTAECJ) 

A\M1)=AV(2)=0, 

AV(3)=1. 

?!  AZTMM=AZIMM/RAO  ? EL E V= EL E V / E AO 

CALL  =CT(AV,AY,PI/2.-£LE\/,l) 

CALL  RCT(AV,AV,AZIMM,3> 

CALL  FOT(AV,AV,-PI/2.+IDcNT(5,J),l) 

DECL-ASIN(AY(3)) 

85  A5-’NC  = ATAN2(Atf(l),AV(2)> 

AS-NC=°I-IOENT (6,J ) ♦GHA-A?'NC 
A 5 "NC-4 M OC ( ASE NC , 2 . *F I ) 

ASCENC=A5ENC*EAO  S DE C LC =D£CL *R A O 
OrLTAE (J)=DELTAEtJ)/RAr 
90  20  FORMAT  (3F10.5) 

AV (1) =AV (2) =0, 

A \l  (3)  =1. 

CALL  RCT(AV,AV,PI/2.-ELEV,l) 

CA  l.L_ROT  ^AV  , A Vj  PI  / 2 •+AZIMM-QECTA£(J)  ,3) 

°5  A SC  = A T AN  2 CAW  ( 1 1 , ■*  A \M  2 ) ) 

CELTAE«J)=A5C*FA0 
TEMPI ( 1) = ASCENC 
TEMPI  (?) =OECLC 
TEMPI ( 3) =DELTAE(J) 

100  _ ASC£NC  = A$CENC/RAO 

OE  CLC=OECLC/RAD ~ 

OELTAE(J)=OELTAE(J)/RAO 
GH(J) =GHA 

IF  ( J. EQ.2. AND. GH( 1) ,EQ.GH(2> > GO  TO  69 
10c  IF(J.EO.l)  GO  TO  69 

OLCH=GH(l)-GH(2)  

GH 1 2) *GH ( 1) 

GHA=GH( 1 ) 

ASCENC-ASCENC-DLGH 
110  TEMPI  ft l = A5CENC*RAD 

69  ACENI J) =ASCENC 
OEO(J)=OcCLC 
rf 

C CALCULATE  THE  AZIMUTH  AND  ELEVATION  OF  THE  COF  OF  THE  CAMERA. 
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11?  r 

av  f i»  = :. o 
av (?)  =:.c 
a 7 ( 3 ? = i . 3 

CALL  RCT(AV,AV,-PI/2.hDcCLC,2> 

12"  PAH  ROT < AV  ,AV,GHA-ASCENC-I0ENT(6, Jl ,3) 

CALI  °0T <AV,AV,PI/2.-I0ENT(5,Jl ,2) 

A.- L ( J ) = A S I N (A  V <31  ) 

AA  L = n L ( J ) *° A C 
A7IM(J)=ATAN2(AV(2),-AV(1>) 

1 2 c AA7IM=A7IM( J) *RAQ 

W=TT£(t,7C)  (TEMPI (I),I=1,3),AAZIM, AAEL.GHA 
’0  ( /,A0X,  ♦RIGHT  ASCENSION  CF  CAMERA  =*  F 09^  3_,  / , 

1-0X,*C  CLINATION  OF  CAMERA  =*F06.3,/, 

?A"y,*QOTATICN  OF  THE  FILM  PLANE  =*F08.3,/, 

133  R///,T2-, ‘AZIMUTH  = *F ft . 3 , T 5- , *t L E V A T ION  = *F" . 3 , T 9- * S I 0- RE AL  TIM 

5F6.3) 
r 

r RE  A3  THE  HEADER  CARD. 

1 7 c =EAC(J,111)  ( TITLE (I ,J) ,1  = 1,8) 

111  FORMAT  ( 3 A 1 0 ) 

WPITE (• ,121)  (TITLMI.J) ,1=1,8) 

121  FORMAT  </,30X,8AlG) 

K=  1 

1 A 3 " 

r AD  TUE  X AND  Y COORDINATES  OF  THE  POINT  RELEASES  IN  THE 

r FILM  PLANE,  AND  CALCULATE  THE  LINE  OF  SIGHT. 

L L = 3 

IAS  125  REA0(J,131)  IX I(L ) , YI( L) , Ir (L ) ,L=1, A)  

LL=LL+i 

171  FORMAT (IX, A ( 2F  8 • A , IA  I ) 

IRtEOF  IJII  180,135 
175  DO  171  L = 1 , A 

15"  IF  (XI  ('  ) .EQ.O. O.ANO.YI (L)  .EQ.O  .0 ) GO  TO  171 

XI  (L)  =-X  ML) 

XX(K|:  XI  L ) 

YY  <K) =vi OL) 

CALL  ANGLE(FL<J)  , XI (L) ,YI (L ) , SIGMA, ETA) 

15?  call  LINSGT(SIGMA,ETA,ASCENC,DECLC,OELTAE(J) ,gha, 

1C  < 1 , K ) ,C(Z,K> , C (3  , K ) ) 

r 

C TH i Sr  STATEMENTS  Sr£  USED  ONLY  HHEN’A  COMPLETE  PRINTOUT  ! IS  “OTSIREO" 

r 

160  GO  TO  55 

Atf(l)=COSIETA)*SIN(SIGMA) 

AV (2) =SINCETA) * S I N (SIGMA, 

A V ( 3) =DOS (SIGMA) 

CALL  ROT ( AV ,AV,-OELTAE ( J)  ,3) 

165  CALL  POT (AV,AV,-PI/2.»OECLC,2) 

CALL  =CT(AV,AV,GHA-ASCENC-I02NT(6,J>,3> 

CALL  RCT(AV,AV,PI/2.-IDENT(5,J),2) 

AA-L=ASIN(AV(3) ) *RAO 

AA7IM=flTAN2(AV(2) ,-AV( 1)) *R  AO 

17  C S I GM  A = 5 1 GM  A * PAD 

eta=fta*rad 
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w 


1 


1?5 


18  0 


195 


WRITE  (5,161)  ID(L)jXICL)  , YI ( L ) ,SIGMA,ETA,AAZIM,AAEL,(C(I),I=1,3> 

161  FOPMAT(5X,I5,3F10.5,F12.5 ,2X,2F7.2,3F10. 5,3X,F7. 3.5X.F7. 3) 

55  K = KU 
171  CONTINUE 
GO  TO  125 
190  N(J)=K-1 

WRIT£(6,80)  N ( J) 

80  FORMAT  (Z,T54,I4,*  NPTS*) 

K=N( J) 

WRITE  (IUNIT, 60)  (XX(JJ),YY(JJ),(C(II,JJ),II=1,3), 

1JJ=1,K> 

60  FORMAT  (2F8.4.3F12.8) 

J_=  J*  1 __  _ 

ENOFILE  IUNIT 
REWINC  IUNIT 


150 


IF<J.£Q.3>  GO  TO  200 
GOTO  6- 
20P  CONTINUE 

WRITE (6,10) 


155 


WRITE  (6,40) 
60  FO°MAT (1HT) 
JSTA= J-l 
KKSUM=0 
SUMX=0. 
SU“X2=C. 


S Y 0 = S Y02  = 0 . 

KPLUS=0 

KMIN=0 


200  C 


WRITE  THE  HEADING  OF  THE  DATA  FILE S,  AND  REWIND  THE  FILES. 


20c 


DO  51  1 = 1 , JST  A 
IFU.EQ.2)  GO  TO  51 
DO  41  J = 1 , JST  A 
IF  ( J.  EO.  I)  GO  TO  >,1 
WRITE (20,225)  SEC 
225  FO^MAMIX^Z.O) 


210 


215 


WR  ITE  (f , 121 > (TITLE(L.I) ,L=1,8) 
WRITE  (6,121)  (TITLEU,  J)  ,1=1,8) 
WRITE (20,111)  (TITLE (L , I) ,L=1, 8) 
WRITE  (20, 111)  (TITLE  (L  , J) ,L  = 1 ,8) 
WRITE  (3,111)(TITLE (L,J) ,L  =1 , 8 ) 

IUNIT  = 10«-I 

JUNlf  = 10  * J 
REWINO  IUNIT 
REWINC  JUNIT 
NN=N( J) 

N N T = N ( I ) 


220 


225 


RE  AO  (JUNIT,  60)  (XX  (JJ)  ,YY  (JJ)  , (C  m,  JJ  ) , 11  = 1 ,3),  JJ  = 1 ,NN) 
CALL  VCTOR (RHO(l , I) ,RHO (1, J) ,0,2) 

CALL  SCALCF (0,0, ZZ ,2) 

Z Z = 1 . /SORT ( ZZ) 

CALL  SCALCF (0, ou, ZZ, 1) 

IF  7 = 0 

J T T = 1 

I N I T = 0 
JT  = JT  T 
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0 AS  SEEN  FROM  ANOTHER  SITE. 

r* 

STCOM='O0 

23* 00  299  KK  = 1 . NNI  

RE  AO ( I UNI T ,60)  X, Y,ICC aiV,II  = l,3» 

U7  00  220  L=l,3 
TE  MP1 (L) =CC  CL ) 

220  CONTINUE 

2A  0 0 

C CALCULATE  THE  INTERSECTION  OF  TWO  LINES  OF  SIGHT. 

p 

ISw  = 0 
TC  = 0 

2<*5  ISW1=C 

* ICCM=3 

FMT  N= 1 . E + 7 
FCCMP=2.E-8 
ZCONP=. 002 
2E  0 Z M I N= 10  0 • 

JT  = JTSAV 
IFJJT.LT. 1»  JT=1 

221  CALL  SCALCFCO,C(l.JT»,TEMF2(l) ,2) 

CTUL  SO  ALCF  (OrTEMPH  1 > , TEMP2  ( 2 I » 2) 

2r 5 CALL  SCALCF(TEMPUl)  ,C(1.JTI  ,TEMP2(3»  ,2) 

AA=1.0-TEMP2(3I*TEMP2(3> 

SL < TE Mp2 ( 1) -TEMP2J2I  * TEMP2  ( 3 » ) / A A 
SL  = I=- CTEMP2<  2I-TEMP2( 1)*TEMP2 (3) > /AA 

CAlCULATE  ERROR  VECTOR- 

CALL  SCALCF<C(1.JT) ,V£MP2<i» ,SLR,1> 

CALL  VICTOR (0,TEMP2(1> ,T€MP2(1) ,2) 

CALL  SC ALCF (TEMPI ( 1) ,TEMP3(1) , SLR 1,1) 

CALL  VCTQR(TEHP2(lt ,TCHP3(1> , TEMP3(1) ,1) 

~C5LL  SCALCF (TEMPT < 11 ,TEMP3(TI , AA,2) 

CALL  V CTOR(OU,TEMPl,AV,3) 

CALL  V-CTOR (OU,C( 1 ,JTI  ,AV2,3) 

CALL  SCALCF(AX,A\/,XXX,2) 

CALL  SCALCF (AV2.AV2, YYY,2> 

XXX  = 1 ,/SQPT  (XXX> 

YYVr'l  ,/SCPT  (YYY» 

CALL  SO  ALCF (AVtAVtXXX.l) 

CALL  SCAlCF (AV2.AV2, YYY,1> 

CALL  SCALCF (AV.AV2 .F ,2) 

F=1.-ABS(F» 

ZZ=SQRT ( A A) 

rP(ISH.FO.l)  GO  TO  260 
IC=IC*t 

IF (ISW1.EO. 01  GO  TO  228 

IF (F.LE. FCOMP.ANO.ZZ.lt. ZMIN>  GO  TO  228 
GO  TO  229 
2 2*  JT  MIN  = JT^ 

FM  TN  = F 
ZMTN=ZZ 


26  5 


27  0 


275 


2?  " 


285 


28 


ISH=1  

229  ir(ICOv.GT.STCOM)  GO  TO  24C 

IF(FMIN.LE.FCOWP.ANO.ZMIN.LE.ZCOMF)  GO  TO  240 
IF  ( IC.GT. ICOM1  GO  TO  231 
232  JT=JT»1 

TF<JT.GT.NN1  GO  TO  240 

GC  TO  ?21  _ 

235  FORMAT)*  NO  SATISFACTORY  MATCH  AFTER* 1 3*  ATTEMPTS*) 
231  CONTINUE 

239  ICCM=TCCM+1 
GO  TO  232 

240  JT=JTMIN 

I S H = 1 

IF (ZMIN.LT.O. 005. ANO.FMIN.LT. FC0MP*5. » GO  TO  221 
HR  ITE (c  * 230  > IC 

GO  TO  299 

26"1  CALL  F>-POS(CC,SLRI,RHO,ACEN,DEC,DELTAE, 

1 GH,FL,XV,YV) 

STCCH=30 

YMIN=SQRT<(XV-XX(JT) )** Z* ( YV- Y Y ( JT)  ) ** 2 ) 

ZZ=ZZ*1000. 

f=f*ic:oooo. 

SY^SYj  + YV-YY  C JT) 

SYC2=SYD2*YMIN 

<~  NEXT  CARO  FOR  SITE  2 FILM  PLANE  PUn 

HR  ITE (3,  265)  XV.YV.KK 
265  F0PHAT(1X,2E11.4,I4) 

JT*AV= JT 
JTT=ZZ 

CALL  SCALCF  (C(1,JT) .TEMPI, SLR, 1) 

CALL  V7CTOR(RHC<l,J> , TEMPI, AV,1) 

ZZ  = 0 • 5 

CALL  SCALCF  (TEMP3 ( 1) ,TEMP3(1)  ,ZZ,1> 

CALL  V'CTOR (AV.TEMP3 (1 ) .TEMP2, 1) 


CONVERT  TO  LAT,  LONG,  ANO  ALT. 


AA  = ASIN(TEMP2<  3) /SORT t TEMF2 (1) **2+TEMP2 (2) ♦* 2+ TEMP2 { 3 ) **2 ) ) 

88  = ATAN2  (TEMP2 (2) , TEMP2<1 )) 

CALL  CCNVPT (2, TEMPI ( 1)  , TEMPI ( 2 I , TEMPI ( 3 I , A A , 86 , TE MP2 ( 1 ) ) 

TEMPI (U  = TEMPI (1) *RAD 
TEMPI (21 = TEMPI (2) *RAO 

WRITE 290)  KK,I, J.YHIN, (TEMP1JII) ,11*1,3) , SLR, SLRI ,F,IC, JT.JTT 
290  F:0';MAT  (9X,  14,  lX,2I3,lX,F8.f  ,1X  ,Fi2.  5,F18.3,Fi7. 3,5X,F10.h,  4X,FiO  . 4 

1,1X,F1C.3,I3,I5,I5) 

WRITE (2  0, 32  0)  (TEMPI (II),  1 1=1, 3) 

320  FORMAT (IX, 3F10 .5) 

K*SUM=KKSUM+1 
IF  (KKSUM. EO.l)  GO  TO  297 
VA-  = TfMPl  ( D-ZLAST 
IF(VAC.LE.O.)  KMIN=KMIN+1 
IF  ( VAP.GT.O .)  KPLUS  = KPLUS*1 
SUMX=S'JMX,  VAR 
SUMX2="UMX?*VAP*»2 
297  ZLAST=TCMP1 (3) 

299  CONTINUE 
7 995  ENOF I t r 20 


29 


141  CONTINUE 
51  CONTINUE 
*AF  FS=*KSUW 

AB-SUNX/FS 

ASIG=SORT (KKSUM*3UMX2-SUMX*  SUMX) /FS 
SYC2  = ISY02/FS) 

SYTrSYT/FS 

350  PRINT  305,A8,ASIG,KKSUM,KNIN,KFLUS,SY02,SYQ 

WRITE (4,10> 

305  FOFMATflHO, ’AVERAGE  ALTITUDE  INCREMENT  = *FA.h/ 
1 » SIGMA  = *F11.6/*  TOTAL  MATCHES  = *16/ 

2 * NC.  ALTITUDE  DECREASES  = *15/ 

3SC  3 * NO.  ALTITUDE  INCREASES  - *16/ 

4 * ABS  Y“  DEVIATION  SITE  2 = *F12.7/' 

5 * AV  Y DEVIATION  SITE2  = *F12.7> 

END  FILE  3 
300  CONTINUE 
3ft!  310  CONTINUE 
STOP 

End 


i 
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SUIROUUNP  ORIENT(V,ASCENC,0-CLC,r:LTAE) 

this  SUBROUTINE  CALCULATES  THU  C=I:NTATICN  OP  THE  CAMERA 

PR  CM  THE  POSITION  OP  A FAIR  OF  3TALS.  ( ALL  ANGLES  Aria  IN  RADIANS) 

DIMENSION  V (3, -) , VC(7, 2) , B( 3) 

IFLAG  = '. 

PI=3. 1-1E02ES 


CALCULATE  DIFFERENCES, 

DC  11  T=l,2 

CALL  V.CTOR(V(l,I)  illll,2»II  » V 0 ( 1 , I ) , 2 ) 

11  CONTINUE 

CALCULATE  AND  SCALE  8. 

CALL  V: CTOR (VD (1, 1 > , VO (1, 2) ,3, 3) 

CALL  SCALCF  (3 ,9, C, 2) 

C = SOP  T (C  ) 

CALL  SCALCF<3,P,1. /C,l) 

CALCULATE  ROTATION  ANGLE. 

lc  CALL  SCALCF  (V(l,l ) ,V(1 .3)  ,C,2> 

CALL  SCALCF  IV (1,1 ) ,3,0,2) 

0=0*0 

COMEGA= (C-O)/ (1.0-0) 

SO-'FGA  = SORT  ( 1 . C -CO  MEGA  * *2  ) 

OMrGA  =AT  AN? (SOMcGA.CCMEGA) 


CALCULATE  ASCENC.CECLC  ANO  OELTAl. 

T1=ATAN2 ( B ( 1 * » E ( 2 ) ) 

T2=SIN(OMEGA/2.0) / CO S ( CMEGA / 2 . C ) 

T3  = ATA>!(B(3)*T2) 

AS0ENC=T3-T1 

0ELTAE=T3*TI 

T1=C0MLGA«-(1.-C0MPGA)*R(3)’B(  3) 
T2=SQRT(1.-T1*T1) 

OECLC  =ATAN2(T1,T2) 

IF  (IFLAG.EC.OI  GO  TO  2E 
RETURN 

25  CALL  V CTORIV (1,3) ,3*VDfl,2) ,3) 

CALL  SCALCF (VO (1, 2), VD( 1,1), SIN (OMEGA) ,1) 
CALL  SO ALCr ( V ( 1 » 3 ) ,3,C,2) 

C = C*(1. -COS  (OMEGA) ) 

CALL  SCALCF(B,V0(1,2),C,1) 

CALL  VrCTOR(VO(l,l),VO(l,2) ,VD(1,1) ,1) 
CALL  SCALCF (V( 1,3) , V0( 1,2 ) ,COS (OMEGA) ,1) 
CALL  V- C TOR  (VO (1,1 ) ,VD(1,2>  ,V0  (1,1) ,1) 

DO  31  J=1  , 3 

IF <ABS( V0(J,1)-V( J,1 )) .GT.G. J50)  GO  TO  5 
31  CONTINUE 
RETURN 

c 00  21  1=1,3 
3 ( T ) = -n ( I ) 

21  CONTINUE 
I FL  AG  = 1 
GO  TC  IE 
END 
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'OP 


SUPRCUUNE  CONVRT  I K , GDL A T ,GDL 0 N , H , GCL A T , GCL 0 N , RHO) 

THIS  SUBROUTINE  CONVERTS  GEODETIC  COORDINATES  TO  GEOCENTRIC 
COORDINATES  WHEN  K=l.  (ALL  ANGLES  ARE  IN  RADIANS) 

WH  N «-?,  THE  INVERSE  TRANSFORMATION  OCCURS. 


DIMENSION  R HO ( 3 1 

CLARK!?  SLLIPSOIO  OF  1866  IS  USED. 

RA  = 6’7S.2C6<, 

RD=635o.5838  

RRt(«.'P/ra>  **2 
EC=1 • -CR 

SELECT  OPTION. 

GO  TO  r ,15)  ,K  

CONVE  PT  GEODETIC  TO  GEOCENTRIC  COORDINATES. 
c GCLCN=GDLON 

0 = CI  A/(cORT(l.  0-EC*  (SIN (GOLAT) > **2) > 

GCLAT  = ATAN(  (RR»H/C)/n.»H/C>»TAN(GDLAT)  ) 

PHO(l ) = (C*H)*COS(GDLAT ) *COS(GOLON) 
=H0(2)=-(C*H>*C0S(GDLAT>*SIN(GCL0N) 

RHO (3)= (C*RM-H)*SIN(GDLAT) 

PETURN 

CONVEPT  GEOCENTRIC  TC  GEODETIC  COORDINATES. 

15  IF(GCLON.LT.D.O)  GCLON=-GCLON 
GDL  CN=GCL  ON 

?-  C=  5 A/ SCR T ( F HO ( 1 1 * RHO ( 1 ) *R HO ( 2 ) *RHO ( 2) ♦RHO (3 ) *PHO(  3 ) ) 

T = AT  AN  (SORT  (RR) * (T  AN(GCLAT)  +C * EC* T A N ( GCL A T > / SORT ( 1 . - EC* < SI N < GC L AT ) 

1>*» ? ) ♦ . 5 * C*  C*  E C*  EC  IN (2. *GCL AT)/ (1 .-EC*SIN (GCLAT ) **  2 > ) ) 

C=TA/C 

H=S0RT(C*C“2. *PA*C*( COS (GCLAT )* COS ( T) ♦ SORT ( RR ) *SI N ( GCL AT > * SI N ( T > ) ♦ 
10  1RA*RA* II . -EC* (SIN( T) **2> ) ) 

GDLAT  = AT  AN(  <C  * S I N ( GCL  A T ) - RA  *SQ  R T (RR)*SIN(T I ) /( C*COS ( GCLAT ) - 
1PA  *COS (T I ) ) 

P E TU°  N _ 

END 
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1 

SUBROUTINE  ANGLE(FL,X, Y, SIGMA, ETA) 

r 

r 

ThIs  subroutine  calculates  the  ANGLES  WHICH 

DEFINE  A POINT 

c 

r 

C 

IN  THE  FILM  PLANE.  ( ANGLES  ARE  EXPRESSEO  IN 

D=SORT (X*X*Y*Y) 

SIGMA=ATAN«0/FL> 

RADIANS) 

ETA=ATAN2(Y»X) 

5 RETURN 

10 

ENO 

1 

SUBROUTINE  LINSGT  < SICMA.ETA, ASCENC,OECLC,OE 

LTAE  , TAU,X,Y,Z  > 

C 

c, 

PERFORM  THE  REQUIREO  ROTATIONS  TO  DETERMINE 

THt_  LINE  OF  SIGHT 

c 

DIMENSION  v/  < 3 » 

PT  = 3. 1“159265358979 
V C 31  = SIN  < SIGMA) 

V ( 1» = COS (ETA) *V<3) 
Y(2I  = SIN(ETA)  »</(3> 

10 

c 

V ( 3) =COS (SIGMA ) 

f 

Q 

POTATE  A SOU  T THE  Z AXIS. 

15 

"AN^-DELTSE 

CALL  ROT(V,V,ANG,3) 

r 

n 

ROTATE  ABOUT  Y AXIS. 

c 

ANG=-PI/2.0+DECLC 

20 

CALL  ROT ( V, V, ANG, 2) 

c 

ROTATE  ABOUT  Z AXIS. 

25 

ANGrTAU-ASCENC 
CALL  ROT ( V , V, ANG , 3 > 

X = V(1 ) 
Y=  Y ( 2 ) 

7 = V ( 3 > 
RETURN 

30 

ENC 
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1 3UER0UTINF  SCALCF (X1,X2,S,N0P) 

r 

C THIS  SUBROUTINE  MULTIPLIES  A VECTOR  XI  BY  4 SCALAR  S WHEN  NLP- 1 , 

r ANH  PLACES  THE  RESULTS  IN  X?. 

5 C IF  NOPr  2 » THE  DOT  PRCOUCT  OF  VECTORS  XI  AND  X2  IS  STCREO  IN  S. 

C 

0 1 NENSI ON  X 1 ( 1 ) ♦ X2  < H 
IF(N0P.GT.2)WRITF(6,5»N0P 
5 F0°MAT(1X,*  AT  SCALCF,  NOP=*.I10) 

10  GO  TO (1 ,3) , NOP 

1 00  2 J=l,3 

X? ( J> = S*X1 ( J> 

_? CONTINUE  _ _ 

RE  TU°  N 

IS  x S=C. 

00  4 J=l,3 
S=S*X1(J)*X2(J» 

4 CONTINUE 

RETURN 

2 C ' END 


1 


r 

5 r. 
c 
0 


10  1 


IF  4 


20 

t- 


SU -TROUT  I NE  VECTORl  XI  ,X  2 , X 3 , NOP  ) 

THTS  SUBROUTINE  PERFORMS  V. CTOR  AOOITION  IF  NCF=1, VECTOR 
SUBTRACTION  IF  NOP=2,  AND  CALCULATES  THE  VECTOR  PRODUCT  IF 
N0F=3.  XI  AND  X2  ARE  THE  INPUT  VECTOPS  AND  X3  IS  THl  RESULTANT 
VECTOR. 

DIMENSION  Xllll ,X2 (1) X3(lt 

GO  T0!1,4,5),NCP 

S=l. 

DO  3 J-1,3 
X 3 t J I =X 1 ( J) +S*X2< J> 

CONTINUE  ____  _ 

RE  T(JPN 
S=-l. 

GO  TO  2 
00  6 J=l,3 
K= J*l-3* ( J/3) 

L?  J»2-T*  ( J/2)  __  

X3(JI=X1(KI*X2(LI-X1(L)*X2(K) 

CONTINUE 

RETUPN 

END 
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SU°ROUTI  Nt  P0T<V1*V3,ARG,NAXIS> 


THIS  SUBROUTINE  PERFOFMS  A TWO  DIMENSIONAL  ROTATION  OF  VECTOR  V 2 
BY  AN  ANCLE  ARG(RAOIANS), WHICH  IS  NEGATIVE  FOR  CLOCKWISE 
ROTATIONS  ANO  POSITIVE  FOR  COUNTERCLOCKWISE  ROTATIONS.  V2  IS  TH£ 
ROT  AT  cO  VECTOR.  NAXIS  IS  THE  AXIS  OF  ROTATION  WHICH  IS  1,2, OR 
3 FOR  NOTATIONS  ABOUT  THE  X,Y  OR  Z AXIS  RES PEC T I VEL Y . 

DIMENSION  Vlfll,V2 (3) , V3C1) 

C = COS  (ARC) 

S=SIN(ARG> 

I=NAXIS+l-3*(NAXIS/3> 

IJ  = NAXTS«-2-3MNAXIS/2> 

V2 (NAXIS»=V1(NAXIS) 

V2(II=V1(I»*C+V1(II)*S 
V2  III )=-Vl(I> *S« VI ( II) *C 

V 3 (1I=02<1) 

V 3 ( 2 • = V 2 ( 2 1 

V 3 ( 3)  =>'2  (3)  _ 

RETURN 

END 


SU3ROUTINE  JULOAY ( JOAY .YEAR, MONTH, DAY) 

JULIAN  0 A Y FOR  MON TH-D A Y- YE  A R IN  20TH  CENTURY 

WHERE  YEAR  = 00-99 
MONTH  = 1-12 

0 A Y = 1-31  __  

DIMENSION  MOSUMI121 
10  INTEGER  YEAR, MONTH  ,OAY , JDAY 

DATA  MOSUM/O,  31,59,90, 120  ,151, 181, 212, 243, 2Z3,  334, 33*./ 

r 

LEAP  = (YEAR  - 1)_  / U _ _ 

JDAY  = 2415019  ♦ 3 6 5 * V E A R * LEAP  ♦ MOSUM (MONTH ) + DAY 
lc  IF  (MCD(YEAR,4) ) 1,2,1 

? IF  ( Y F AR | 3,1,3 

t IF  I MONTH  - 2)  1,1,4 

4 JOAY  = JOAY  ♦ 1 

1 RETURN 

20  E NO 
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SUBROUTINE  FPPQS  f CC ♦ SLRI  .RHO . ACEN  .DEC . DELTAE  . 

1 GH,FL,XV,YV» 

DIMENSION  CC(3>,RHQ(3,3> ,ACEN(3) ,0EC(3) ,0ELTA£13J  . GH(3), 
i Fim 

DIMENSION  TEMPI  ( 3 1 *T£MP2L3»  , T£MPJ  (11 

ic=o 

IK=1 

PI=3. 141 59265358979 

1=1  

J=2 

CALL  SCALCF(CC, TEMPI, SLRI, 1) 

CALL  VICTOR (TEMPI, RHO C 1,1), TEMPI, 1) 

AA  = ASIN( TEMPI  (3)/SQRT(TEHFim**2*TEMPl(2»**2*TEMPl(3>**2>  ) 

6B=ATAN2 (TEMPI (2» .TEMPI (1)) 

CALL  CONVRT(2,GLA,GLO, AT. ,AA , Be, TEMPI > 

AC  1 = 0 . I 

AC2=0.C 

00  1121  KI=1, 3 
AC1=RH0(KI,J)**2*AC1 


AC2=TEMP1 (K 11**2* AC2 

1121  CONTINUE  

CALL  VECTOR (TEMPI, RHO( 1, J), TEMPI, 2) 

R0T=0ELTAE( J) 

CALL  RLNSGT (SIG, ETTA, ACEN (J), DEC (J) ,ROT, 

1GH(J) .TEMPI (11 .TEMPI (2) .TEMPI  I 31) 


00  = FL  (J>*TAN(SIG> 

IF (ETTA. GE, 0.0. AND. ETTA .L E.PI/2. I GO  T 0 1122 

IF (ETTA. GE. PI/2. 0 . AND. ETTA. LE. PI)  GO  TO  1123 

IF  (ETTA.LE.O. O.ANOjETTA.GE. -PI/2.1  GO_ TO  1124 

ETT  A=PI*ETT  A 

TEMP2 (IK)=-DO*COS(ETTA> 

TEMP3 (IK) =-00*SIN ( ETTA ) 

GO  TO  21  

1122  TEMP2(TK»=OO*C0S(ETTA) 

TEMP3(IK)=0D*SIN(ETTA)  

GO  TO  21 
1124  ETTA=-ETTA 


TEMP2(IK»=0D*C0S(ETTA» 

TEMP3(IK)=-00*SIN(ETTA)  

GO  TO  21 

1123  ETTA  = PI-ETT A 

TEMP2 (IK ) =-00* COS ( ETTA  > 

TEMP3(IK»=00*SIN(ETTA) 

21  CONTINUE 

YV=TEMP3(1»  

X V = TEM°2 (1) 

RETURN  

ENO 


